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Abstract. Real-time computation of time-dependent quantum mechanical problems are 
presented for nuclear many-body problems. Quantum tunneling in nuclear fusion at low 
energy is described using a time-dependent wave packet. A real-time method of calculating 
strength functions using the time-dependent Schrodinger equation is utilized to properly treat 
the continuum boundary condition. To go beyond the few-body models, we resort to the 
density-functional theory. The nuclear mean-field models are briefly reviewed to illustrate its 
foundation and necessity of state dependence in effective interactions. This state dependence 
is successfully taken into account by the density dependence, leading to the energy density 
functional. Photoabsorption cross sections in 238 U are calculated with the real-time method for 
the time-dependent density-functional theory. 



1. Time-dependent approaches to quantum problems in nuclear physics 

In this paper, we will report several applications of the real-time calculations. Although we 
concentrate our discussion on the nuclear physics problems, the real-time approaches are useful 
in many fields of many-body quantum systems. Before discussing applications, let us start 
from giving some reasons why we adopt the real-time calculations, instead of a more prevalent 
time-independent method. 

1.1. Stationary solutions 

It is customary to solve a quantum mechanical problem in a form of stationary equations, such 
as the eigenvalue problem of the Schrodinger equation: 



If a direct solution of this problem is achieved, the solution provides a wave function of each 
energy eigenstate. In principle, we may calculate all the physical observables using these wave 
functions. Amount of the computational task is roughly proportional to number of eigenstates 
|$ n ) to be calculated. Thus, for instance, when we are interested in the excitation spectra and/or 
strength functions in a wide range of energy, it should be computationally very demanding. 
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1.2. Time- dependent solutions 

In contrast, a real-time solution of the time-dependent Schrodinger equation is, in general, a 
linear combination of energy eigenstates \& n )' 

\^(t))=J2 C ne- lE ^ n ). (2) 

n 

The number of eigenstates and their coefficients c n involved in Eq. ([2]) are simply determined 
by the initial state. A trivial but important feature of the time-dependent solution, Eq. ([2]), is 
that the state consists of many eigenstates and distribution of their weights, c n , can be chosen 
as we wish, by selecting the initial wave packet. Therefore, we may obtain physical quantities 
for states at a variety of energies from a single time evolution of the wave packet, using the 
Fourier decomposition into different energies. Since the computational task does not strongly 
depend on the weights c n , the real-time calculation for the time-dependent equation could be 
computationally more efficient than solving the time-independent eigenvalue equation ([T]) to 
obtain many eigenstates. 

1.3. Boundary condition for continuum states 

Another practical advantage of the real-time calculation is the fact that the method does not 
require the boundary condition for the description of the unbound continuum states. The 
solution of Eq. ([1]) for an unbound state at positive energy E n cannot be uniquely determined 
without the asymptotic boundary condition. In contrast, the solution of the time-dependent 
equation, Eq. ([2]), is unique for a given initial state. In this sense, the boundary condition is 
naturally given by the time-dependent solution. 

To illustrate this fact, let us consider scattering of two particles A and B, which are initially 
apart from each other, outside of the range of interaction V . These particles may be composite. 
The Hamiltonian Hq describes two free particles, A and B, at the initial channel. The total 
Hamiltonian of the system is given by H = Hq + V. The scattering solution for the energy 
eigenstate can be formally written as [T] 

l * <+,(E » > = ( 1 + ¥3^")l*< E »- < 3 > 

where the initial (incident) state is an eigenstate of Hq, Hq\<&(E)) = E\&(E)). The Green's 
function, (E — H + irj)^ 1 , contains the infinitesimal ir\ that determines the scattering boundary 
conditions. The state of Eq. ([3]) is obtained by acting the M0ller's wave operator on the initial 
state. 

= 17(0,-00)1$^)) = (l -/ dt'UMe-^Vjit')) ME)), (4) 

where U(t,t') = e iH t e -iH(t-t') e -iH t' and _ e iH ty e -iH t Herej the initial state |$(£)) is 
given by a plane wave for the relative motion between A and B. Then, we need the convergence 
(adiabatic switch-off) factor e" 7 ''*'', in order to remove the interaction between A and B at the 
limit of t — > — oo. This trick is necessary when we treat the energy eigenstate, because the 
particles A and B are overlapped in the plain-wave solution and interacting each other. In 
contrast, the infinitesimal factor r\ becomes unnecessary [T] if we use a wave packet, |$ pac ), 
instead of the eigenstate \$>(E)) in Eq. (JU. Since the wave packet is spatially localized, we can 
construct the initial state in which the wave packets of A and B are far apart, not interacting each 
other (Fl^pac) = 0). The scattering wave packet, which is a superposition of many eigenstates, 
can be calculated by propagating the initial wave packet. We may stop the time evolution at 
a finite period of time, because the interaction takes place only in the finite time. Once the 



wave packet is scattered away outside of the interacting region, the state is governed by the 
Hamiltonian Hq. After the scattering, the proper scattering asymptotic behavior should be 
automatically imposed for the wave packet. 

1.4- Intuitive picture of quantum dynamics 

The time-dependent description has not only these advantages in practical computation, but 
also provides an intuitive picture of the quantum dynamics. It is very helpful for our insight 
into complex systems to visualize movement of the wave packet. This will be demonstrated in 
the following applications. 

1.5. Nuclear many-body problems 

The nucleus is a self-bound quantum system which presents a rich variety of phenomena. It 
is composed of fermions of spin 1/2 and isospin 1/2, called nucleons (protons and neutrons), 
interacting with each other through a complex interaction with a short-range repulsive core [2]. 
Thus, the nucleus is a strongly correlated system. Furthermore, it is highly "quantum", which 
could be quantified by the magnitude of the zero-point kinetic energy relative to the potential 
energy, (T)/(V) ~ 1. This is roughly unity in nuclear systems: (T)/(V) ~ 1, which indicates 
that the nucleus has a very strong quantum nature. Even for the liquid helium, this ratio is 
smaller than the nuclear case. 

Remarkable experimental progress in production and study of exotic nuclei requires us to 
construct a theoretical model with higher accuracy and reliability. Extensive studies have been 
made in the past, to introduce models and effective interactions to describe a variety of nuclear 
phenomena and to understand basic nuclear dynamics behind them [21 [3]. Simultaneously, 
significant efforts have been made in the microscopic foundation of those models. For light 
nuclei, the "first-principles" large-scale computation, starting from the bare nucleon-nucleon 
(two-body & three-body) forces, is becoming a current trend in theoretical nuclear physics. It 
is computationally very challenging, because, as we mentioned above, the nucleus is a highly 
quantum, strongly correlated Fermionic many-body system, interacting via the complex and 
"singular" nuclear force. Although these ab-initio-type approaches are still limited to nuclei 
with very small mass number, they have recently shown a significant progress. These issues are 
addressed by other contributions [H El [6] • 

In this paper, we show a few nuclear-physics problems: The low-energy nuclear fusion 
reactions, the nuclear response in the continuum, and the time-dependent density-functional 
theory in the linear regime. Here, we would like to emphasize again that, although a particular 
model/theory is adopted for each calculation, the concept of the real-time method is quite 
general and applicable to other models and to other subfields of physics as well. 

2. Real-time calculation of sub-barrier fusion reaction 

In this section, we present the time-dependent wave-packet method and discuss tunneling 
dynamics in sub-barrier fusion process El [91 [THl CD]. To illustrate the essential idea, we 
start from a simple two-body problem. 

2.1. Two-body time- dependent wave-packet model 

A simple potential model of the fusion reaction is constructed as follows: The Hamiltonian 
for the relative motion between two nuclei is given by (h = 1) 
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This leads to the radial Schrodinger equation in each partial wave: 



where ip(R,t) = ^2lm( u l(R, t) /r)Yi,M (R) • The real potential, V(R), consists of the repulsive 
Coulomb and attractive nuclear potentials. When the approaching two nuclei pass through the 
Coulomb barrier, this is regarded as the fusion. To calculate its probability, it is convenient to 
use the imaginary potential, iW(R), that is non-zero only inside the Coulomb barrier. Then, we 
define the fusion by the flux loss caused by iW(R). We assume here the collision of 10 Be [Z = 4) 
on 208 Pb (Z = 82). W(R) and the nuclear part of V(R) are assumed to be in the Woods-Saxon 
form. 

The most common way to solve this problem is to find a stationary solution with a fixed 
energy E by integrating the time-independent radial Schrodinger equation. Then, in the 
asymptotic region, one can compare the flux of incoming and outgoing Coulomb waves, ul{R) ~ 

C { £u ( ~ ] {R) + Cl ut u [ + ] {R), to obtain the fusion probability, P L = (|Q n | 2 - \Cf*\ 2 )/\Cf\ 2 . 
Results obtained in this way are shown by circles in the bottom panel of Fig [2] for L = 0. 

We demonstrate an alternative approach: the time-dependent wave-packet method. The 
initial wave function is of the Gaussian form 

u L (R, t = 0) = e -iKR e -j 2 (R~Ro) 2 ^ ( 7 ) 

where Rq indicate the initial position of the Gaussian center. The parameters, K and 7, are 
arbitrary, but should be chosen so as to cover the interested energy region. In the top panel of 
Fig. [JJ the solid line shows the initial wave packet for L = 0, the head-on collision with average 
relative momentum K. In the present calculation, we set K 2 /2fi = 28 MeV and calculate 
the time evolution up to tf = 10 MeV -1 ~ 2000 fm/c. Time evolution of the wave packet is 
calculated by recursive operation of the small-time evolution operator, U(At) = e~ l ^ tH , that 
is simply approximated with the fourth-order Taylor expansion, U(At) ~ Xm=o( — iAtH) n /n\. 
The wave packet collides with the Coulomb barrier in the second and third panels in Fig. [TJ 
Then, the wave packet is partially reflected back. In the bottom panel, the final wave packet 
contains only outgoing waves. The missing part of the wave packet passes through the barrier, 
entering the inner region and disappearing because of the absorbing imaginary potential, iW(R). 
One can see, in the second and third panels, that a part of the wave packet actually penetrates 
the Coulomb barrier. This incoming component vanishes in the fourth and fifth panels. 

Now, let us show the energy profile of the initial and final wave packets. We calculate the 
energy distribution of the wave packet at time t. 

f L (E;t) = (u L (t)\5(E-H)\u L (t)). (8) 

This quantity can be also calculated using the time propagation technique [7J. The dashed line 
in the top panel of Fig. [2] indicates the calculated energy distribution of the initial wave packet, 
fL=o(E,ti). This is in a Gaussian form whose centroid is about 39 MeV that corresponds to 
the given kinetic energy of 28 MeV. The difference of 11 MeV is the Coulomb potential energy 
at Rq ~ 40 fm. The final distribution, fo(E,tt) shown in the solid line in Fig. (top panel), 
is calculated from the wave packet at the bottom of Fig. [TJ We have fo(E,tf) ~ for E > 42 
MeV. These high energy parts end up the fusion. On the other hand, the flux remains almost 
invariant (fo(E,tf) ~ fo(E,ti)) for E < 36 MeV, which suggests that the fusion cross section 
is negligible in this low energy region. These arguments can be quantified by calculating the 
fusion probability, 

Prfm _ h(E,ti)-f L (E,t f ) 
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Figure 1. Snap shots of time-dependent wave 
functions, Re ul=o{v). The panels are ordered 
in time from the top to the bottom. The 
dashed (dotted) line indicates a shape of the real 
(imaginary) part of the potential. See text for 
details. 



Figure 2. Top panel: Energy distribution 
of the initial (dashed) and final wave 
packets (solid line) corresponding to the top 
and bottom panels in Fig. [TJ respectively. 
Bottom: Calculated fusion probability 
(L = 0) as a function of energy. The 
solid line is calculated with the time- 
dependent wave-packet method, while the 
circles indicate results obtained from time- 
independent (fixed energy) solutions of the 
Schrodinger equation. 



This is shown in the bottom panel of Fig. [2] by solid line. The result calculated with the 
time-dependent wave-packet method perfectly agrees with that of the fixed-energy calculation 
(circles). This means that the time-dependent method provides an alternative method of 
accurate calculation for fusion cross section. 

2.2. Three-body time- dependent wave-packet model 

Let us upgrade the model to three bodies and discuss effects of a valence neutron on the 
fusion. There are many theoretical and experimental investigations on this, but the conclusion 
is somewhat elusive yet. Especially, the effect of weakly-bound valence neutrons with a spatially 
extended wave function, which is often called "neutron halo" , is controversial. There were many 
arguments that the weakly-bound neutron enhances the fusion cross section, especially at sub- 
barrier energies |12 | 113 1 IT^ l H5| , I16j. However, as is shown in the followings, we have reached the 
opposite conclusion [71 151 191 H01 ITT] . 

We assume that the projectile is a bound state of a neutron (n) and a core nucleus (C), 



and the target (T) is treated as a point particle. In this work, we neglect the intrinsic spin of 
neutron. The time-dependent Schrodinger equation of the three-body scattering is given as 

ift|^(R, r, t) = {-^l ~ + Kc(r) + V ct {Rot) + KtKt)} *(R, r, t), (10) 

where we denote the relative n-C coordinate as r and the relative P-T coordinate as R. The 
reduced masses of n-C and P-T motions are m and fi, respectively. The n-C potential, V n c(r), 
is assumed to be real with the Woods-Saxon form. This potential must produce a bound state 
of the neutron around the core. The depth of the potential V n c is varied according to the 
orbital energy simulating a projectile nucleus with a tightly-bound to weakly-bound neutron. 
The core-target potential, Vct(-Rct), contains an imaginary part iW(Rct) that simulates the 
fusion. The n-T potential, VnT^TVr), is taken to be real. The imaginary potential is present 
only between the core and target nuclei. When the core and target become close enough, the 
wave function will vanish, to be counted as fusion. Therefore, since the final destination of 
the neutron is irrelevant, the present calculation does not distinguish complete and incomplete 
fusion. The total fusion probability (sum of complete and incomplete ones) is calculated below. 

The wave-packet is constructed using the partial wave expansion in the body-fixed frame 
|17l I18j. This calculation is equivalent to the one in the space-fixed frame. However, the body- 
fixed frame has an advantage for calculations of states with non-zero total angular momentum 
J 7^ (J = L + l). In the body-fixed frame, channels are characterized by the magnetic quantum 
number (the projection of J on the z-axis) and the angular momentum I conjugate to the 
angle 9 between r and R. Coupling between different Q channels is present only for AO = ±1, 
caused by the Coriolis term. In the present paper, we focus our discussion on the head-on 
collision (J = 0). See Ref. [10] for the complete calculation. The J = state is characterized 
by Q = and L = I, given by 

^°(R, r, 6; t) - 'g r - ' > ^ fltcos t). (1!) 
1=0 r 7r 

Here, we need to truncate the space by setting the maximum value of the partial waves I. To 
calculate the Coulomb breakup process, we may take Z max = 4. However, when the n-T potential 
is included, nuclear breakup and the neutron transfer may take place, then, much larger Z max is 
necessary to obtain convergent results; Z max ~ 70 [H[T0]. Since the present calculation is simply 
based on the discretization of the radial coordinate, it is straightforward to include these high 
partial waves. This is an advantage over the continuum-discretized coupled-channel calculations 

In this model, the projectile is assumed to have a simple structure that a neutron sits on the 
2s orbital in the n-C potential V Q c- <fto( r ) = A/l/4vryo( r )/ r - The initial wave packet at time 
t = is prepared in the same manner as we have shown in Sect. 12.11 Namely, the P-T relative 
wave function is given by the boosted Gaussian wave packet whose centroid is at Rq = 25 fm. 

Ul (R,r,t = 0) = e-^ R e-^ R - R ^y Q {r)5 m . (12) 



2.2.1. Case (1) Well-bound projectile: First, we investigate fusion reaction of projectile with 
a well-bound neutron. For this case, we found that the fusion cross section is identical to the one 
without the valence neutron. Thus, the valence neutron essentially sticks to the core nucleus, 
and they move together as a single nucleus. This may be naturally expected. 

However, there are exceptions, when the condition of energy matching is satisfied, which 
means that energies of the neutron orbital in the projectile and target are approximately 
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Figure 3. Fusion probability (J = 0) as a function of incident energy for halo Be (e 



-0.6 

MeV) on 209 Bi. Dotted line indicates the probability for the projectile without the valence 
neutron. Dashed and solid lines indicate fusion probabilities calculated with and without V n T 
potentials, respectively. 




Figure 4. Time-dependent density distribution of the valence neutron, p(r, 9; t) in the collision 
process of 11 Be on 209 Bi. The panel at the left (right) end shows the initial (final) wave packet. 
The center of each panel, r = 0, indicates the position of the core, and the right is the direction 
for the target, 9 = 0. 



degenerate, we observe a substantial enhancement of sub-barrier fusion cross section. The fusion 
probability is strongly correlated with the neutron transfer probability and sensitive to V^t- This 
may be understood in terms of adiabatic dynamics of the valence neutron. See discussion in 
Ref. 0. 

2.2.2. Case (2) Weakly-bound projectile (Neutron halo): Now, let us discuss fusion of weakly 
bound n Be on 209 Bi, namely effects of the neutron halo. As is mentioned above, there are 
theoretical arguments that the weakly-bound neutron enhances the sub-barrier fusion cross 
section |12l [131 [14"1 [15| [16] . We studied the same problem in a fully microscopic framework with 
the real-time method. 

Calculated fusion probabilities at the head-on collision ( J = 0) are shown in Fig. [3l The 
fusion probability is slightly suppressed by the presence of the halo neutron. It seems that the 
fusion dynamics for the neutron-halo projectile are mainly determined by the Coulomb breakup 
process. Time evolution of the neutron density distribution during the collision process is shown 
in Fig. HI We plot the following quantity in the x-y plane with x = r cos 9 and y = r sin 9 for 



the head-on collision of J = 0; 



p(r,9;t) = J dR\t J=0 (R,r,6;t)\ 2 . (13) 

The direction of the target always corresponds to 9 = 0, which means that the target nucleus 
approaches from the right side on the x-axis, then returns back to the right. Figured suggests the 
followings: When the core nucleus is decelerated by the target Coulomb field, the halo neutron 
behaves like a spectator, keeping its incident velocity and leaving the core nucleus. This yields 
the Coulomb breakup of the projectile, then reduces the collision energy between the core and 
the target. This spectator picture of the fusion dynamics for halo nuclei nicely explains why the 
total fusion cross section does not increase even though the radius of the neutron-halo nucleus 
is so large [19]. 

Our conclusion is that the presence of the halo neutron suppresses the fusion cross section, 
irrelevant to the collision energy. This contradicts other former predictions [12|. \T3[ [T4"l \T5[ flB] . In 
fact, the model adopted in the coupled-channel calculation in Ref. [16] is very close to ours. We 
have found that the difference comes from the truncation of the model space |10| . Namely, the 
result in Ref. [16J did not reach the convergence with respect to the number of the partial wave 
l max in Eq. (jlip . If we set l max < 4 as is done in Ref. |16j . we also obtain an enhanced fusion 
cross section at sub-barrier energies. However, the cross section decreases as £ max increases, and 
finally it becomes even smaller than the two-body result without the valence neutron. 



3. Real-time calculation of strength functions in the continuum 

In this section, we discuss a method of calculating the strength distribution, in particular, the 
El strength distribution, from the time-dependent wave function. 
The El strength distribution in the continuum is defined by 



dB(El-E) 



dE 



m,m' 



JdE'5(E-E') (<S>( + \E';m')\M(El;m)\<$> ) 2 ■ (14) 



Here, M {El; m) is the electric dipole operator with the magnetic quantum number m. The final 
states \&( + \E; m)) should have the proper outgoing asymptotic form when unbound channels 
are open. The construction of the final state in the continuum is a difficult task in general, 
especially for many-body non-spherical systems. An alternative way of avoiding this difficulty is 
the time-dependent description. Let us show a simple example as an illustration of the method. 



3.1. Illustrative example 

The El strength function of Eq. (|14p can be written in a time-dependent form: 
dB(El; E) l T „\-^.« f°° jj-ME+irit/x. i * ,U m . ^\ „-itH 



dE 
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This is easily obtained by the fact that the time integration simply produces the Green's function 
(E — H + ir/) _1 whose imaginary part is nothing but the delta function 5(E — H). Therefore, 
we can calculate the El strength distribution as follows: (1) Construct the initial state by 
operating the El operator to the ground state, \^(t = 0)) = M(El;m)\$o)- (2) Calculate its 
time evolution and its overlap with the initial state. (3) Its Fourier transform leads 

to dB (El; E) / dE . In this time-dependent formulation, we do not need to construct the final 
states with the proper boundary condition. This is automatically taken into account in the 
time evolution of wave packets. When the El operator excites the ground state into unbound 



continuum states, the wave packet will decays in time and the integrand of Eq. ()15p will 

vanish at t — > oo. Thus, we do not need the convergence factor 

Let us discuss application of the time-dependent method to a simple single-particle model. 
We again use the Woods-Saxon potential of V n c in Sect J2,2l to describe the n Be nucleus as a 
two-body system of a 10 Be core and a neutron. The neutron is assumed to sit on the 2s orbital. 
Multiplying it by the El operator with the recoil charge, we construct the initial state which is 
shown in the top panel of Fig. We solve the time-dependent Schrodinger equation similar to 
Eq. ^ with L = 1. The imaginary potential iW(R) is active at R > 50 fm to absorb outgoing 
waves. The time evolution is shown in the following panels in Fig. [SJ The state with L = 1 
propagates into r — > oo in time. After a certain period, the state no longer has an overlap 
matrix element with the initial state. Then, we stop the calculation. From this time evolution, 
we obtain El strength distribution (dotted line) shown in Fig. There is a strong low-energy 
peak known as the threshold effect. This peak is drastically reduced if we change the energy 
of the 2s orbital. For instance, the solid line in Fig. [6] shows the result for the case that the 
energy of the 2s orbital is —2 MeV. Experimental data [20|, [21j [22] indicate a similar shape but 
somewhat smaller strength. 

This method can be extended into many-body problems in a straightforward manner. 
However, its computational task drastically increases as the number of particles increases. Thus, 
we resort to the density functional model of nuclei in the next section. 

4. Density-functional approach to nuclei 

The density functional theory is a leading theory for describing nuclear properties of heavy nuclei 
and perhaps the only theory capable of describing all nuclei and nuclear matter with a single 
universal energy density functional. In the nuclear physics, it is often called the self-consistent 
mean-field model, because of a historical development based on the Brueckner-Hartree-Fock 
theory and introduction of the effective interaction. However, it should be noted that there are 
both fundamental and practical differences between the naive mean-field theory and the mean- 
field model of the nucleus which is, in fact, conceptually analogous to the density-functional 
theory. 

I briefly review basic properties of nuclei and show that those properties cannot be understood 
by a simple mean-field model. Here, the saturation property plays a key role. 

4-1- Nuclear saturation in the mean-field model 

Nuclei are known to be well characterized by the saturation property. Namely, they have an 
approximately constant density po ~ 0.17 fm~ 3 , and a constant binding energy per particle 
B I A ~ 16 MeV0 In this section, I show that the nuclear saturation property has a great impact 
on nuclear models. Especially, it is inconsistent with the independent-particle model of nuclei 
with a "naive" average (mean-field) potential. 

There are many evidences for the fact that the mean-free path of nucleons is larger than 
the size of nucleus. In fact, the mean free path depends on the nucleon's energy, and becomes 
larger for lower energy [2]. Therefore, it is natural to assume that the nucleus can be primarily 
approximated by the independent-particle model with an average one-body potential. The 
crudest approximation is the degenerate Fermi gas of the same number of protons and neutrons 
(Z = N = A/2). The observed saturation density of po ~ 0.17 fm" 3 gives the Fermi 
momentum, kp ~ 1.36 fm -1 , that leads to the Fermi energy (the maximum kinetic energy), 
Tp = k F /2M w 40 MeV. 

1 This is the extrapolated value for the infinite nuclear matter without the surface and the Coulomb energy. 
The observed values for finite nuclei are B/A ~ 8 MeV. 
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Figure 5. Time evolution of a state (real 
part of the wave function) initially excited 
by the dipole operator. The vertical scale 
for r = 10 ~ 100 fm is magnified by a factor 
of two compared to that for r = ~ 10 fm. 
The dashed line indicate the shape of the 
p-wave potential. 
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Figure 6. Calculated El strength 
distribution for 11 Be. See text for 
explanation. 



First, I show that the independent-particle model with a constant attractive potential V < 
cannot describe the nuclear saturation property. It follows from the simple arguments. The 
constancy of B/A means that it is approximately equal to the separation energy of nucleons, S. 
In the independent-particle model, it is estimated as 

Sf* B/A rj ~{T F + V). (16) 

Since the binding energy is B/A ~ 16 MeV, the potential V is about —55 MeV. It should be 
noted that the relatively small separation energy is the consequence of the significant cancellation 
between kinetic and potential energies. The total (binding) energy is given by 

where we assume that the average potential results from a two-body interaction. The two kinds 
of expressions for B/A, Eqs. (|16p and (|17p . lead to Tp ~ — 5V/4 ~ 70 MeV, which is different 




from the previously estimated value (~ 40 MeV). Moreover, it contradicts the fact that the 
nucleus is bound (Tp < |V|). 

To reconcile the independent-particle motion with the saturation property of the nucleus, 
the nuclear average potential should be state dependent. Allowing the potential V% depend on 
the state i, the potential V should be replaced by that for the highest occupied orbital Vp in 
Eq. (|16p . and by its average value (V) in the right-hand side of Eq. (|17p . Then, we obtain the 
following relation: 

V F « (V) + T F /5 + B/A. (18) 

Therefore, the potential Vp is shallower than its average value. 

Weisskopf suggested the momentum-dependent potential V, which can be expressed in terms 
of an effective mass m* [23] : 

Vi = U Q + uA. (19) 

Actually, if the mean-field potential is non-local, it can be expressed by the momentum 
dependence. Equation (|19p leads to the effective mass, m* /m = (1 + Ux/Tf)' 1 . Using Eqs. 
(fT7j|) . (fT8"j) . and (HHP, we obtain the effective mass as 

m* f3 5B11" 1 , 

+ o-T7^ -0.4. (20) 
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Quantitatively, this value disagrees with the experimental data. The empirical values of the 
effective mass vary according to the energy of nucleons, 0.7 < m*/m < 1, however, they are 
almost twice larger than the value in Eq. (|20p . As far as we use a normal two-body interaction, 
this discrepancy should be present in the mean-field calculation with any interaction, because 
Eq. (|20p is valid in general for a saturated self-bound system. Therefore, the conventional 
models cannot simultaneously reproduce the most basic properties of nuclei; the binding energy 
and the single-particle property. 

This problem can be solved in the density functional theory. In nuclear physics, basically the 
same solution was interpreted in terms of the density-dependent interaction. The variation of 
the total energy density functional E[p] with respect to the density contains "re-arrangement 
potential", dV e s[p]/dp, which appear due to the density dependence of the effective force V^ffp]. 
These terms turn out to be crucial to obtain the saturation condition. Now, the expression for 
the total energy, Eq. (|17j) . should be modified to include the re- arrangement effect. This resolves 
the previous issue, then provides a consistent independent-particle description for the nuclear 
saturation. This is important to understand basic nuclear properties and to describe nuclei in a 
universal framework. 

Intensive studies in nuclear density functional models in recent years have produced numerous 
results and new insights into nuclear structure \24\ 125] . However, it is impossible to review all 
of them in this paper. Thus, I will present results of our study with the (time-dependent) 
density functional approach, on the shape phase transition in the rare-earth nuclei and on the 
photoabsorption cross section in 238 U. 



4-2. Density-functional theory for superfluid nuclei 

A modern energy functional for nuclei is a functional of many kinds of density, such as kinetic 
r(r) and spin-orbit density J(r). In addition, it is known that the heavy nuclei in the open-shell 
configurations show characteristic features of the superfluid systems [3J. To describe superfluid 
properties produced by pairing correlation among nucleons, we need to add the pair (abnormal) 
density K,(r). These densities are collectively denoted as p in the followings. Variation of the total 



energy, E[p], leads to the following equation, known as the Hartree-Fock-Bogoliubov equation 
in nuclear physics [3]: 

ff[* J **]|* M )=^|* M ), (21) 

where and are quasi-particle energies and 
states, respectively [3]. is composed of two 

components; the upper \U^) and the lower one |V^). 
The solution of Eq. (|21|) defines the normal density 
p(f) = Eii^W^M) the pair density K(f) = 
t^u(^)^u(0) an d other densities at the ground 
state. Since /if'S/, \&*] depends on these densities, 
Eq. (f2Tj) must be solved in a self-consistent way. 
Minimization of the energy density functional may 
lead to a spontaneous breaking of symmetry. An 
example is given in Fig. [7] for Nd isotopes [26]. 
The intrinsic quadrupole moment calculated with 
the Skyrme functional of SkM* is compared with 
the experimental data. At N = 82, the nucleus 
at the ground state is spherical Qq = 0, while for 
N = 86 ~ 92, the deformation gradually develops. 
The observed ground-state deformations deduced 
from the transition probability B(E2; 2 + — > + ) are 

nicely reproduced. Note that there are no adjustable parameters in this calculation. 
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Figure 7. Calculated (open squares) and 
experimental (filled) intrinsic quadrupole 
moment for Nd isotopes 



4-3. Time- dependent density-functional theory for superfluid nuclei 

For a description of the time-dependent phenomenon, we need to extend the energy functional 
to include the time-odd densities, such as the spin density s(r, t) and the current density j(r, t). 
Now, all the densities are time dependent. The time-dependent version of Eq. (|2ip is formally 
written as 

i^ t \* fl (t)) = HMt),y*(t)}\*M- (22) 

This is known as the time-dependent Hartree-Fock-Bogoliubov equation in nuclear physics. 

The computation of Eq. (|22|) is very demanding because we need to calculate the time 
evolution of all the quasi-particle states The number of them is same as the dimension 

of the single-particle model space which is much larger than the particle number. Although we 
proposed a feasible approach to its linear regime, known as the finite amplitude method (FAM) 
J2Z1 [23 [291 GSB EI], the FAM is based on a time-independent formulation. In this section, we 
present another approximate treatment of Eq. (j22[) to provide time-dependent equations. 

The approximate equations alternative to Eq. (|22p is called the canonical-basis TDHFB |32j. 
Assuming the diagonal property of the pair potential in the canonical pair of states k and k, 
Eq. (|22p can be approximated by the following set of equations. 

i J^k(t)> = (h(t) - Vk (t))\Mt)), * J^*(*)> = (Kt) ~ V- k (t))\h(t)), (23a) 
ij t Pk(t) = K k (t)A%(t) - K* k (t)A k (t), (23b) 



These basic equations determine the time evolution of the canonical states, and \(() k (t)}, 

their occupation, p k (t), and pair probabilities, K k (t). The real functions of time, ijk(t) and ^(t), 



are arbitrary and associated with the gauge degrees of freedom. The time-dependent pairing 
gaps, Afc(i), are defined in the same manner as the BCS pairing gap [3] except for the fact that 
the canonical pair of states are no longer related to each other by time reversal. It should be 
noted that the quantities in the Cb-TDHFB equations, (p, k, A), are not matrixes but are only 
their diagonal elements with a single index for the canonical states k. 

The Cb-TDHFB equations, (|23p . are invariant with respect to the gauge transformation with 
arbitrary real functions, 6 k (t) and k (t). 

\<Pk)^e w ^\<p k ) and |^)^ e %«|<^) (24) 
Kk ^ e -i(W)+e- k (t)) Kk and Afe e -m(t)+e- k (t)) Ak (25) 

simultaneously with 

Vk(t) ^ m (t) + and m(t)^Vk(t) + ^-- 

It is now clear that the arbitrary real functions, r] k (t) and rj k (t), control time evolution of the 
phases of \(j> k (t)), |0k(*)), n k (t), and A k (t). 

In addition to the gauge invariance, the 
Cb-TDHFB equations possess the following 
properties. 

(i) Conservation law 

(a) Conservation of orthonormal prop- 
erty of the canonical states 

(b) Conservation of average particle 
number 

(c) Conservation of average total energy 

(ii) The stationary solution corresponds to 
the HF+BCS solution. 

(iii) Small-amplitude limit 

(a) The Nambu-Goldstone modes are 
zero-energy normal-mode solutions. 

(b) If the ground state is in the normal 
phase, the equations are identical 
to the particle-hole, particle-particle, 
and hole-hole RPA with the BCS 
approximation. 

For numerical calculations, we extended 
the computational program of the TDHF in 
the three-dimensional (3D) coordinate-space 
representation [33] to include the BCS-type pairing correlations. The ground state is first 
constructed by the HF+BCS calculation. Then, we add a weak impulse electric dipole field 
Vei(*) to the ground state. VEi(t) is chosen as 

V E1 (t) = -r,Ff 1 6(t), Ff 1 = Yl (Ne/A)n- £ (Ze/A)n, (26) 

i:protons i:neutrons 

where i = (x, y, z). We solve the Cb-TDHFB equations in real time and real space. To obtain the 
El strength function, we calculate the time evolution of the expectation value (t)\F E1 \ 1 <i' (t)) 
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Figure 8. Calculated (solid line) and 
experimental (symbols) photoabsorption [33] 
cross sections for 238 U. The dashed and dotted 
lines indicate the contribution from the dipole 
oscillation parallel and perpendicular to the 
symmetry axis, respectively. 



under the external field Vgi(i) with small rj. Then, the strength function S(E1;E) can be 
obtained by the Fourier transform of (^>(t)\F E1 \^(t)) with an artificial damping factor with 
r = 1 MeV: {^(t)\F E1 \y(t)) -)• {^(t)\F E1 \^(t)) e - rt / 2 . 

The calculated El strength distribution is transformed into the photoabsorption cross section 
and shown in Fig. [8l for the 238 U nucleus. The ground state of 238 U is deformed in the prolate 
shape with the axial symmetry. Thus, the photoabsorption peak is split into two peaks: the 
lower peak represents an oscillation along the symmetry axis (K = 0) and the higher one 
corresponds to that perpendicular the the symmetry axis (K = 1). The total cross section is 
the sum of these, which well agree with the experimental data. This can be another indication 
of the deformation of the 238 U nucleus. 

The computational task for solving Eq. (|23p is significantly smaller than that for solving Eq. 
([22]) . This is due to the fact that the number of the calculated canonical states k, k is roughly 
the same order as the particle number (N ~ 10 2 ), which is much smaller than the dimension of 
the model space (M ~ 10 4 - 10 5 ). 

5. Summary 

The real-time calculation is a useful tool to calculate many-body dynamics, especially when we 
are interested in bulk properties in a wide range of energy, and when the complicated continuum 
boundary condition is necessary to construct the energy eigenstates. We have demonstrated its 
capability and usefulness in several examples, including the fusion reaction and the response in 
nuclei. The methodology is quite universal and applicable to a variety of sub-fields in many-body 
quantum physics. 
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